#!/usr/bin/perl -w
#bins sequences by %GC in 3rd codon position and counts occurrences of each codon within each bin
#output format: (tab separated)
#first line: number of sequences in each bin
#second line: average %GC3 of sequences in each bin
#64 following lines: codon followed by number of times codon is counted in each sequence bin

$fasta = $ARGV[0];#fasta formatted DNA file (sequence all on 1 line)
open(FASTA,$fasta);

#$output = "CIGC3-" . $fasta;
$output = $ARGV[1];
open(OUTPUT,">$output");

$binsize = 5;#non-overlapping bins: each bin will later be combined with next bin to form partially over-lapping bins
$bins = 100/$binsize;#number of bins

@SeqCount = ();#sequences in each bin
until (scalar @SeqCount == $bins){push(@SeqCount,"0");}

@SumGC3 = ();#sum %GC3 for averaging
until (scalar @SumGC3 == $bins){push(@SumGC3,"0");}

@Codons = ();#counts for each codon
until (scalar @Codons == 64){push(@Codons,"0");}

@CodonBins = ();#counts for each codon in each %GC3 bin
$count = 0;
until ($count == $bins){push @CodonBins, [@Codons];++$count;}

while ($sequence = <FASTA>){
	unless ($sequence =~ ">"){#collect sequences, ignore header lines
	
		chomp $sequence;
		$length = length $sequence;
		$codons = $length/3;
		
		if ($codons == int($codons)){#if complete reading frame
			$GC3 = 0;
			$CdnEnd = 2;#start with last base of first codon
			until ($CdnEnd > $length){#until last codon: check 3rd pos
				$base3 = substr($sequence,$CdnEnd,1);#last base of current codon
				if ($base3 eq "G" or $base3 eq "C"){++$GC3;}#if last base is GC
				$CdnEnd = $CdnEnd + 3;
			}#until last codon
			
			$perGC3 = 100*($GC3/$codons);#%GC in 3rd codon pos
			$GCbin = int($perGC3/$binsize);#GC bin number (array position)
			
			++$SeqCount[$GCbin];#count sequences in each bin
			$SumGC3[$GCbin] = $SumGC3[$GCbin] + $perGC3;#sum %GC3 values in each bin for averaging
			
			#count occurrences of each codon in sequence:
			$CdnStart = 0;
			until ($CdnStart + 3 > $length){#until last codon
				$codon = substr($sequence,$CdnStart,3);
				
				if ($codon eq "TTT"){++$CodonBins[$GCbin][0];}
				elsif ($codon eq "TTC"){++$CodonBins[$GCbin][1];}
				elsif ($codon eq "TTA"){++$CodonBins[$GCbin][2];}
				elsif ($codon eq "TTG"){++$CodonBins[$GCbin][3];}
				elsif ($codon eq "CTT"){++$CodonBins[$GCbin][4];}
				elsif ($codon eq "CTC"){++$CodonBins[$GCbin][5];}
				elsif ($codon eq "CTA"){++$CodonBins[$GCbin][6];}
				elsif ($codon eq "CTG"){++$CodonBins[$GCbin][7];}
				elsif ($codon eq "ATT"){++$CodonBins[$GCbin][8];}
				elsif ($codon eq "ATC"){++$CodonBins[$GCbin][9];}
				elsif ($codon eq "ATA"){++$CodonBins[$GCbin][10];}
				elsif ($codon eq "ATG"){++$CodonBins[$GCbin][11];}
				elsif ($codon eq "GTT"){++$CodonBins[$GCbin][12];}
				elsif ($codon eq "GTC"){++$CodonBins[$GCbin][13];}
				elsif ($codon eq "GTA"){++$CodonBins[$GCbin][14];}
				elsif ($codon eq "GTG"){++$CodonBins[$GCbin][15];}
				elsif ($codon eq "TCT"){++$CodonBins[$GCbin][16];}
				elsif ($codon eq "TCC"){++$CodonBins[$GCbin][17];}
				elsif ($codon eq "TCA"){++$CodonBins[$GCbin][18];}
				elsif ($codon eq "TCG"){++$CodonBins[$GCbin][19];}
				elsif ($codon eq "AGT"){++$CodonBins[$GCbin][20];}
				elsif ($codon eq "AGC"){++$CodonBins[$GCbin][21];}
				elsif ($codon eq "CCT"){++$CodonBins[$GCbin][22];}
				elsif ($codon eq "CCC"){++$CodonBins[$GCbin][23];}
				elsif ($codon eq "CCA"){++$CodonBins[$GCbin][24];}
				elsif ($codon eq "CCG"){++$CodonBins[$GCbin][25];}
				elsif ($codon eq "ACT"){++$CodonBins[$GCbin][26];}
				elsif ($codon eq "ACC"){++$CodonBins[$GCbin][27];}
				elsif ($codon eq "ACA"){++$CodonBins[$GCbin][28];}
				elsif ($codon eq "ACG"){++$CodonBins[$GCbin][29];}
				elsif ($codon eq "GCT"){++$CodonBins[$GCbin][30];}
				elsif ($codon eq "GCC"){++$CodonBins[$GCbin][31];}
				elsif ($codon eq "GCA"){++$CodonBins[$GCbin][32];}
				elsif ($codon eq "GCG"){++$CodonBins[$GCbin][33];}
				elsif ($codon eq "TAT"){++$CodonBins[$GCbin][34];}
				elsif ($codon eq "TAC"){++$CodonBins[$GCbin][35];}
				elsif ($codon eq "TAA"){++$CodonBins[$GCbin][36];}
				elsif ($codon eq "TAG"){++$CodonBins[$GCbin][37];}
				elsif ($codon eq "CAT"){++$CodonBins[$GCbin][38];}
				elsif ($codon eq "CAC"){++$CodonBins[$GCbin][39];}
				elsif ($codon eq "CAA"){++$CodonBins[$GCbin][40];}
				elsif ($codon eq "CAG"){++$CodonBins[$GCbin][41];}
				elsif ($codon eq "AAT"){++$CodonBins[$GCbin][42];}
				elsif ($codon eq "AAC"){++$CodonBins[$GCbin][43];}
				elsif ($codon eq "AAA"){++$CodonBins[$GCbin][44];}
				elsif ($codon eq "AAG"){++$CodonBins[$GCbin][45];}
				elsif ($codon eq "GAT"){++$CodonBins[$GCbin][46];}
				elsif ($codon eq "GAC"){++$CodonBins[$GCbin][47];}
				elsif ($codon eq "GAA"){++$CodonBins[$GCbin][48];}
				elsif ($codon eq "GAG"){++$CodonBins[$GCbin][49];}
				elsif ($codon eq "TGT"){++$CodonBins[$GCbin][50];}
				elsif ($codon eq "TGC"){++$CodonBins[$GCbin][51];}
				elsif ($codon eq "TGA"){++$CodonBins[$GCbin][52];}
				elsif ($codon eq "TGG"){++$CodonBins[$GCbin][53];}
				elsif ($codon eq "CGT"){++$CodonBins[$GCbin][54];}
				elsif ($codon eq "CGC"){++$CodonBins[$GCbin][55];}
				elsif ($codon eq "CGA"){++$CodonBins[$GCbin][56];}
				elsif ($codon eq "CGG"){++$CodonBins[$GCbin][57];}
				elsif ($codon eq "AGA"){++$CodonBins[$GCbin][58];}
				elsif ($codon eq "AGG"){++$CodonBins[$GCbin][59];}
				elsif ($codon eq "GGT"){++$CodonBins[$GCbin][60];}
				elsif ($codon eq "GGC"){++$CodonBins[$GCbin][61];}
				elsif ($codon eq "GGA"){++$CodonBins[$GCbin][62];}
				elsif ($codon eq "GGG"){++$CodonBins[$GCbin][63];}

				$CdnStart = $CdnStart + 3;
			}#until last codon
			
		}#if complete reading frame
	}#unless header: seq line
}#while reading fasta file
close FASTA;

@AverageGC3 = ();#average %GC3 for each %GC bin
#combine each bin (except last) with next bin to form partially overlapping bins:
#0-10,5-15,10-20, etc, up to 90-100
$bin = 0;
until ($bin == $bins - 1){
	$Seq1 = $SeqCount[$bin];#number of seqs in small non-overlapping bin
	$Seq2 = $SeqCount[$bin + 1];
	$Seqs = $Seq1 + $Seq2;#number of sequences in larger bin
	$Sum1 = $SumGC3[$bin];#sum of %GC3 in small non-overlapping bin
	$Sum2 = $SumGC3[$bin + 1];
	if ($Seqs > 0){$AveGC3 = ($Sum1 + $Sum2)/$Seqs;}else {$AveGC3 = "NA";}#take average %GC3 for bin unless there are no sequences in bin
	push(@AverageGC3,$AveGC3);
	print OUTPUT $Seqs,"\t";#print number of seqs in each bin
	++$bin;
}

print OUTPUT "\n";

foreach $Average (@AverageGC3){print OUTPUT $Average,"\t";}#print average %GC3 for each bin
print OUTPUT "\n";
@CodonNames = ("TTT","TTC","TTA","TTG","CTT","CTC","CTA","CTG","ATT","ATC","ATA","ATG","GTT","GTC","GTA","GTG","TCT","TCC","TCA","TCG","AGT","AGC","CCT","CCC","CCA","CCG","ACT","ACC","ACA","ACG","GCT","GCC","GCA","GCG","TAT","TAC","TAA","TAG","CAT","CAC","CAA","CAG","AAT","AAC","AAA","AAG","GAT","GAC","GAA","GAG","TGT","TGC","TGA","TGG","CGT","CGC","CGA","CGG","AGA","AGG","GGT","GGC","GGA","GGG");
$CodonNum = 0;

#foreach codon and each %GC3 sequence bin: print number of codon within bin
until ($CodonNum == 64){
	print OUTPUT $CodonNames[$CodonNum],"\t";
	$curbin = 0;
	until ($curbin == $bins - 1){
		$codon1 = $CodonBins[$curbin][$CodonNum];
		$codon2 = $CodonBins[$curbin + 1][$CodonNum];
		$codons = $codon1 + $codon2;
		print OUTPUT $codons,"\t";
		++$curbin;
	}#until print all bins
	print OUTPUT "\n";
	++$CodonNum;
}


close OUTPUT;
exit;